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ABSTRACT 

A study of the statistics of cosmological black-hole jet sources is applied to EGRET blazar data, 
and predictions are made for GLAST. Black-hole jet sources are modeled as collimated relativistic 
plasma outflows with radiation beamed along the jet axis due to strong Doppler boosting. The 
comoving rate density of blazar flares is assumed to follow a blazar formation rate (BFR) , modeled 
by analytic functions based on astronomical observations and fits to EGRET data. The redshift and 
size distributions of gamma-ray blazars observed with EGRET, separated into BL Lac object (BL) 
and flat spectrum radio quasar (FSRQ) distributions, are fit with monoparametric functions for the 
distributions of the jet Lorentz factor T, comoving directional power l' e , and spectral slope. A BFR 
factor « 10 x greater at z > 1 than at present is found to fit the FSRQ data. A smaller comoving rate 
density and greater luminosity of BL flares at early times compared to the present epoch fits the BL 
data. Based on the EGRET observations, w 1000 blazars consisting of 800 FSRQs and FR2 radio 
galaxies and » 200 BL Lacs and FR1 radio galaxies will be detected with GLAST during the first 
year of the mission. Additional AGN classes, such as hard-spectrum BL Lacs that were mostly missed 
with EGRET, could add more GLAST sources. The FSRQ and BL contributions to the EGRET 
7-ray background at 1 GeV are estimated at the level of « 10 - 15% and « 2% - 4%, respectively. 
EGRET and GLAST sensitivities to blazar flares are considered in the optimal case, and a GLAST 
analysis method for blazar detection is outlined. 

Subject headings: AGNs: blazars — black holes — gamma-ray bursts 



1. INTRODUCTION 

Population studies of black-hole jet sources, which in- 
clude blazars, gamma-ray bursts, and microquasars, are 
difficult because of the unknown emission processes and 
beaming patterns of the relativistic jets. Moreover, the 
density and luminosity evolution of black-hole jet sources 
through cosmic time is uncertain. Here we develop a 
method to treat the statistics of black-hole jet sources 
using the 7-ray data alone. Although the focus of this 
study is radio galaxi es and blazars. the m ethod can also 
be applied to GRBs (|Le fc Dermer|[2006l ). 

The interest in population statistics of blazar sources 
is that an accurate determination of source density 
evolution is needed t o identify parent populations 
(|Urrv fc PadovariiH l995). to chart black-hole formation 
and growth throug h out t he history of the u n iverse 
dBottcher fc Dermerl 120021: ICavaliere fc D'Elial 120021 : 
iMaraschi fc Tavecchioi 12001 . and to assess the contri- 
bution of black-hole jet sources t o the 7-ray background . 
The isotropic 7-ray background (|Sreekumar et al]|1998f) 
consists of an extragalactic 7-ray background (EGRB) 
and an uncertain contribution of quasi-isotropic Galac- 
tic 7 rays produced, for example, by Compton-scattered 
radia tions from cosmic-ray electrons ([Strong et al.ll200"ol 
12001 . 

Soon af ter the recognition o f the 7-ray blazar class with 
EGRET (|Fichtel et a lj 1 19941) ■ 7-rav b lazar population 
studies were undertaken. iChiang et all (|1995f ) performed 
a (V/Vmax) analysis assuming no density evolution and 
showed that luminosity evolution of EGRET blazars was 
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implied by the data. With a larger data set and using 
radio data to ensure the sa mple was unbiased in regard 
to redshift determination, IChiang fc Mukherieej (1998) 
again found that evolution was required. They obtained 
best-fit values through the maximum likelihood method 
that gave an AGN contribution to the EGRB at the level 
o f sa 25%. 

iStecker fc Salamo"nl ([1996) postulated a radio/7-ray 
connection, and tried to correct for the duty cycle and 
7-ray spectral hardening of flaring states. They found 
that essentially 100% of the EGRET EGRB arises from 
unresolved blazars and AGNs. They did not, however, 
fit the blazar redshift distribution to provide a check on 
their model. The crucial underlying assumption of this 
approach, whi ch has been developed in further detail 
in rec ent work (jGiommi et a l. 2006; Nar umoto fc Totanil 
2006), is that there is a close connection between the 
radio and 7-ray properties of blazars. Because a large 
number of E GRET blaza r s (FS RQs) are found in the 5 
GHz, > 1 Jv iKiihr et al.l ((19811) catalog, a radio-/7 ray 
correlation is expected, but is not found in 2.7 and 5 
GHz monitoring of EGRET 7-ray blazars (|Mucke et al.l 
1997). X-ray selected BLs are also not well-sampled in 
radio surveys. Studies based on correlations between the 
radio and 7-ray emissions from blazars may therefore be 
based on a questionable foundation. It is therefore neces- 
sary to distinguish between the very different properties 
and histories of FSRQs and BLs and their separate con- 
tributions to the EGRB. 

A detailed physical model to treat blazar statistics 
that avoids any radio/7-ray blaz ar correlation was de- 
veloped by iMiicke fc Po hl (2000). Blazar spectra were 
calculated assuming an injection electron number index 
of —2. Distributions in injected particle energy into BL 
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and FSRQ jets were considered in the modeling. The 
indices in the injected energy distributions were taken 
from the luminosity functions of FR1 and FR2 radio 
galaxi es which, according to t he blazar unification sce- 
nario (|Urrv fc Padovariil Il995f ). are the parent popula- 
tions of BLs and FSRQs, respectively. A simple descrip- 
tion of density evolution is given in the form of a cutoff at 
some ma ximum redshift z max . Depending on the value 
°f z max, Miicke & PohJ (|2000D concluded that as much 
as w 40 - 80% of the EGRB is produced by unresolved 
AGNs, with w 70 - 90% of the emission from FRls and 
BLs. 

Here we also treat a physical model f or FSRQs and 
BLs. This is similar to the study of iMiicke fc Pohll 
( 2000), though different in a number of important ways. 
No detailed radiation modeling is employed, but rather 
we use mean spectral indices as measured by EGRET. 
Various model forms describe the rate density of blazar 
flares, but only single values of luminosity and T fac- 
tor are considered for each of the FSRQ and BL classes. 
A mono-luminosity function for blazars means that the 
range in apparent powers is kinematic, arising from the 
different, randomly oriented jet directions. This simple 
blazar model is highly constrained when fitting to data, 
even given the freedom to consider different redshift- 
dependent analytic forms for the blazar formation rate 
(BFR). In the case of detailed fits to blazar data, the 
assumption of no luminosity evolution can be and, for 
BLs, is relaxed. 

We use parameter sets that give acceptable agreement 
to the EGRET data on blazar redshift and size distri- 
butions to make predictions that will be tested with 
GLAST, 2 and to estimate the blazar contribution to 
the EGRB. The connection of the properties of 7-ray 
blazars to blazars detected at radio, X-ray, and other 
wavelengths can be used to determine the accuracy of 
models that assume a radio/7-ray connection. 

Section 2 describes the sample of EGRET blazars used. 
The equations for the analysis are presented in Section 
3, and results of the parameter study are described in 
Section 4. Predictions for GLAST and estimates for 
the EGRB from unresolved blazars are presented in Sec- 
tion 5, and we summarize in Section 6. Sensitivities of 
EGRET and GLAST to blazar flares, expressions for op- 
timal sensitivities, and a discussion of a GLAST analysis 
strategy for blazar populations are given in the Appendix 
A. The self-absorption frequency is derived in Appendix 
B. 

2. SAMPLE 

Crucial to making a proper comparison of a model to 
EGRET blazar data is to choose a sample that is un- 
biased with respect to exposure and background. For 
example, exposure to the region around 3C 273 and 3C 
279 was much longer than average over the lifetime of 
EGRET, so blazars found in these pointings would be 
detected to much smaller flux thresholds than on av- 
erage. Likewise, blazars in the vicinity of the galactic 
plane would have to be much brighter than high-latitude 
sources to be detected above the diffuse galactic 7-ray 
emission. 

The 18 month all-sky EGRET survey (|Fichtel et all 
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Fig. 1. — EGRET observations of the redshift distributions of 
blazars, separated into FSRQ (solid) and BL Lac (dotted) popu- 
lations, are shown by the histograms. Smooth solid curve shows 
a monoparametric FSRQ blazar fit with BFR IR,4 (see Fig. 3), 
r = 10, p = 3.4, EC statistics, and l' e = 2.5 X 10 39 ergs s" 1 sr _1 . 
The dotted curve shows the monoparametric blazar fit for the BL 
data for a blazar model with V = A, p = 3.0, l' e (z) = 6 X 10 42 z 1! ' 5 
ergs s — 1 sr — 1 , synchrotron/SSC statistics, and tibl(z) <x z -9 / 4 . 

I1994I ). which ran from 1991 May to 1992 November, 
had roughly uniform exposure over all parts of the sky. 
Thirty-eight AGN identifications were reported in this 
catalog. Additional analysis of the Phase 1 data , as 
reported in the 3EG catalog (|Hartman et al.l ll999D . re- 
vealed numerous additional detections of AGNs during 
Phase 1. 

We have identified all the high-confidence blazars listed 
in the Third EGRET catalog that also appear during 
the 18 month all-sky survey, during which all parts of 
the sky receive roughly uniform exposure. The sample 
we use consists of 60 high-confidence gamma-ray blazars, 
consisting of 14 BLs and 46 FSRQs. We exclude sources 
within 10° of the Gal actic p l ane, and use source catalogs 
(jPadovani fc Giommi|[l99a iPerlmanl fl996l ) to establish 
BL identifications. The integral photon number fluxes 
4>{> E), in units of ph(>100 MeV) cm" 2 s -1 , were used 
to construct the FSRQ and BL size distributions. 

Ta ble 1 lists the sources from the Third EGRET cat- 
alog (Hartma n et al.lll999l ) used in this study and their 
classifications. The 7-ray blazar sample, binned by red- 
shift and grouped into BLs and FSRQs, is plotted in Fig. 
[TJ Fig. [2] shows the FSRQ and BL size distributions. As 
shown in Appendix A, the on-axis EGRET sensitivity for 
a two-week pointing to high-latitude blazars is 1O _8 0_8 
ph(> 100 MeV) cm~ 2 s" 1 , with 0_ 8 ^ 15. This is in 
accord with the dimmest blazars detected with EGRET 
during the two-week pointings. The corresponding vF v 
threshold flux i s f f hr = E 2 cf> s (E), where the source flux 
is given by eq. (jAlljl . Thus 

fl hr = 1.6 x 10" 12 (a ph - 1)0_ 8 ergs cm" 2 s' 1 (1) 

for measurements at E > 100 MeV (e > 200), where a p h 
is the number spectral index. For = 15, f* hr = 2.4 x 
10~ n (a p /j — 1) ergs cm -2 s _1 , noting that a p h — (p+l)/2 
for Thomson, synchrotron, and SSC processes, where p 
is the number index of the assumed power-law electron 
distribution. As shown in Appendix A, an energy range 
that depends on the source spectral hardness should be 
used for optimal detection significance of blazars with 
GLAST or EGRET. 

The BLs and FSRQs display very different 7-ray prop- 
erties. The BLs are less numerous and closer, with an 
average redshift (z) ~ 0.37. For the FSRQs, (z) ~ 1.11, 
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Fig. 2. — EGRET observations of the size distribution of blazars 
in terms of the integral photon fluxes at photon energies E > 100 
MeV versus model size distributions for FSRQ and BL blazars, 
with parameters given in the caption to Fig. 1. Also shown are 
the integral flux sensitivities for EGRET in the pointing mode for 
a two-week observation (over « 1/24 of the full sky) and for 
GLAST in the scanning mode for one year (over the full sky). 
The line labeled "—3/2" has the slope of the size distribution of 
monoluminous sources uniformly distributed in flat space. The 
inset shows the size distribution of EGRET blazars in terms of 
apparent blazar luminosities in the energy range 100 MeV - 5 GeV. 

with a tail on the FSRQ distribution reaching to z ~ 2.3. 
There are « 5x as many FSRQs as BLs per unit peak 
flux in the EGRET range 0„ 8 « 25 - 100. The FSRQs 
are ~ 1 - 2 orders of magnitude more luminous than the 
BLs (inset to Fig. [5]). The apparent powers of the FSRQs 
are as large as « 10 50 ergs s _1 , compared to BLs, which 
are typically < 10 48 ergs s -1 . 

More detailed studies of the redshift and size distribu- 
tion can be made with the EGRET data (R. Romani, 
private communic ation, 2006) based on surveys to iden- 
tify 7-ray blazars (jSowards-Emmerd et al.ir2005h . These 
studies also show that the mean redshift of FSRQs is > 1, 
and that both the BL and FSRQ populations have tails 
to high redshifts. Th e highest redshift 7- ray blazar can- 
didate has z = 5.47 (|Romani et alj |2004). Thus we can 
be sure that the FSRQ BFR extends to at least z > 5. 

The Third EGRET catalog gives the two-week average 
fluxes for blazars during Phase 1. In comparison with 
models, the flaring timescale of the blazars is given by 
a rate density, which relates to the source density by 
a flaring timescale and duty cycle factor. We take this 
flaring timescale to be « 2 weeks for the observer, and 
set the duty cycle equal to unity. If the duty cycle is less 
than unity, then the source space density must be greater, 
leaving the calculation of the intensity from unresolved 
sources independent of these factors. 

3. ANALYSIS 

We employ a simplified version of the standard m odel 
for blazars considered bylDermer fe Gehreld (|1995h (see 
also ISikora et all 119971 : iDermer fc Davis! l2000f )"~A rel- 
ativistically moving plasmoid ejected from a black-hole 
engine has accelerated within it, either through inter- 
nal or external shocks or otherwise, a power-law distri- 
bution of quasi-isotropic ultra-relativistic electrons with 
number index p. In its proper frame, the plasmoid is as- 
sumed to entrain a randomly oriented magnetic field with 
mean strength B. The nonthermal electrons emit syn- 
chrotron or Thomson radiation that is Doppler boosted 



by the effects of the relativistic motion. For the cal- 
culations shown here, we use the synchrotron beam- 
ing factor oc 5^ +p ^ 2 for synchrotron and synchrotron 
sclf-Compton (SSC) processes, where the Doppler fac- 
tor <5d = [r(l — /3cos(9)] -1 , and the external Comp- 
ton (EC) beaming factor oc <5q +p (|Dermerlll995l ). which 
also holds for scattering in the Klein-Nishina regime 
(|Georganopoulos et alj|200ll) . Here T is the bulk Lorentz 
factor of the plasma blob, (3 = VT — T~ 2 , a — (p— 1)/2 = 
oiph — 1 is the energy spectral index of the radiation, and 
6 is the angle between the jet and line-of-sight directions. 
Continuous outflow sc enarios produce a beamin g pattern 
weaker by one power (jLind fc Blandfordlll985[ L but the 
single plasmoid approximation is more applicable to the 
flaring blazars. 

3.1. Cosmology of Nonthermal Sources 

Consider a power-law distribution of electrons with 
low- and high-energy cutoffs, so that the total number 
distribution of nonthermal electrons within the plasmoid 
is described by 



n'm = Ki~ p h (7; 71,72) 



(2) 



Note that p is the number index of the emitting electrons, 
and could be different from the injection index if cooling 
is important. Normalizing to the total comoving electron 
energy W' e — m e c 2 j l d-f 7 N' e (^) implies 
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Coolin g can introduce a b reak in the electron spec- 
trum jMiicke fc Pohll I2OOO0 . but the EGRET spectra 
of bright blazars are well fit by a single power law 
(|Mukheriee et al.|[l997t ). indicating that the single power- 
law approximation is adequate for a treatment of blazar 
statistics. 

The vF v nonthermal synchrotron radiation spectrum 
for a comoving isotropic power-law distribution of elec- 
trons entrained in a randomly oriented magnetic field is 
given in the 5-function approximation by the expression 
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where the luminosity distance for a flat ACDM universe 
is 



d L (z) 



Hr 



(1 + z) / dz 



(5) 



Ub = B 2 /8ir is the magnetic- field energy density in the 
jet plasma, b = B / B cr , and B cr = m 2 c 3 /eh is the critical 
magnetic field. 

The vF v spectrum of jet electrons that Thomson scat- 
ter an external quasi-isotropic monochromatic radia- 
tion field with stationary (explosion)-frame dimension- 



less photon energy e* = 
energy density ET* is 
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(|DermeH (|1995| ); see lDermer fc Schhck eiser 2002, for ex- 
pressions describing the Thomson-scattered accretion- 
disk radiation fields). Restriction to the Thomson regime 
implies that e z < l/(8e«), so that a target 5 eV UV ra- 
diation field would display effects from the onset of the 
KN decline in the cross section at E > 6 GeV/ (1 + z). A 
more accurate treatment for GLAST analyses will have 
to consider the effects of the KN decline on the statistics. 

The vF v synchrotron self-Compton (SSC) radiation 
spectrum in the 5-function approximation is 
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S c = 



where the Compton-synchrotron logarithm 

hx(a m axj 'drain), a m ax = mmffr'Y^ e77i , e ' ~ X ), C l min = 

max(67i,e77|) 1 a nd e' = e z /5 B (jGouldl Il979t 
IDermer et al.lll997D . The SSC process has a similar de- 
pendence as the synchrotron process — though a curva- 
ture in the spectrum is produced by S c — but with a dif- 
ferent coefficient that depends on the physical size of the 
radiating plasma. 

In this formulation, the radiating plasma is spherical in 
the comoving frame, with volume V b ' = 4irr' b 3 /3, where 
the blob radius 
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and t var (s) = 86400t(day) = 10 3 t3 is the measured vari- 
ability time scale. For FSRQs measured with EGRET, 



r' b < 



2.6 x 10 16 (r/10)t(day) 



cm ~ 0.01t(day)(r/10) pc 



(l + z)/2 

and for BLs measured with EGRET 

r' b < 2.4 x 10 14 (r/4)i 3 cm ~ 10 _4 i 3 (r/4) pc . 

For these blob sizes, the synchrotron radiation at GHz 
frequencies is likely to be heavily self-absorbed, as shown 
in Appendix B. This calls further into question the use 
of any radio/7-ray correlation. 
We rewrite eqs. ((4J — (J7J) as 
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where a v — (3 — p)/2 is the vF v spectral index, the 
directional comoving luminosity 
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the beaming factor index 



EC 
(10) 



(p + 5)/2, synchrotron, SSC 



.P + 3, 



EC 



and U Bcr = S 2 ,/8tt. 

The event rate per sr (or the directional event rate) for 
bursting sources in a ACDM cosmology is 



dN _ c f°° d\{z) hco m (z) 

dn ~H Q J z (1 + z )3 y/n m (i + z )3 + n A 



, (11) 



where fi com (z) is the rate density of sources at redshift 
z (see iDermerl I2006L for a detailed d erivation). From 
the WMAP data (jSpergel et all l2003t ). we take n m = 
0.27, n A = 0.73, and Hubble's constant H = 72 km s" 1 
Mpc -1 . The directional event rate (i.e., number count) 
of steady sources is 
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(12) 



and n com (z) is the differential source density. 

In an integral formulation for bursting sources (|Dermerl 
11992ft . the observed directional event rate above the vF v 
spectral flux threshold fl hr of the telescope is 
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Here fP roc is the vF v flux for the process under consid- 
eration, with synchrotron, EC, and SSC fluxes written 
as eq. (|9|), fi — cos 9, and the blazar jet luminosity is 
charaterized by l' e . For a better treatment of detector re- 
sponse, one should calculate a photon-energy integration 
over effective area, rather than describing a 7-ray tele- 
scope by a vF v flux sensitivity ]f lT at a single photon 
energy e. 

For mono-parameter S- function distributions of p, T, 
and I' we have 



dN(> fl hr ) _ 2c 



dn 



= - — / dz 



d\(z) h com (z) 



Ho Jo (1 + zf ^n m {i + zf + n A 



d[i / 4f 6 *[/ e -yr oc (Cp.r,M)], (i4) 
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where the various parameters specifying emission prop- 
erties, for example, W' e , B, r' b , and [/*, are found in the 
directional power l' e . In this expression, we include a 
factor of 2 for a two-sided jet. 

3.2. Peak Flux and Size Distribution 

It is trivial to perform the integration over df e in eq. 
(fT4|) . which places a limit on the allowed values of fi. 
Only values of cosine angle /1 > fi give detectable fluxes, 
where 

fi(z,T,l' e ,p, q ,e,f e ) = f, = I [l-I(gl) 1/<? ] . 

The vF v flux size distribution of blazars per sr per s is 
therefore given by 

dN = 2c r°° ^ d\(z) n com (z)[l-max(-l,A)] 

dn U) Ho Jo (1 + zf y/n m (i + zf + n A 

(16) 

for two-sided blazar jet sources. 
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3.3. Redshift Distribution 

The directional redshift distribution of sources with 
vF v flux / e > fl hr is simply given by 



dN 



thr\ 



2c 



4 



dzdVL 

(z) n com (z) 



[1 — max(— 1, fi)] 



Ho (i + zf ^n m (i + z)i + n A 

The threshold limitation prescribed by jl in eq. (|15|) now 
has / e — > /*'"". The relation between f* hr and the inte- 
gral photon number flux variable </>_g depends on p and 
a v as given by eq. (JTJ) . 

3.4. Parameters 

Eq. (J9j) shows that a 7-ray blazar can be detected from 
redshifts of order unity at the level / e when the param- 
eters satisfy the condition 



r«(2e) Q " 
d\(z = l) 



(18) 



which holds when T 3> 1. If a significant fraction of 
the EGRB originates from blazars, then p = 3.2. For 
synchrotron, SSC, or Thomson processes, this implies 
a photon number index a p h — 2.1, or a v = —0.1. 
EGRET observations of blazar spectral indices show that 
a ph = 2.03 ± 0.09 for BLs and a ph = 2.20 ± 0.05 for FS- 
RQs, with an average sp ectral index of a p h = 2.15 ±0.04 
(jMukheriee et alj 119971 ). The evidence from EGRET 
that a p h is harder for BLs than for FSRQs means that the 
vF v peaks of the 7-ray components of BLs are typically 
at higher energies than for FSRQs. Moreover, observa- 
tions suggest that the flaring state spectra are harder 
than the q uiescent emission, at le ast in the case of PKS 
0528+134 (jMukheriee et al.lll996l ). 

Taking a„ = (p = 3) and noting that d^(z = 1) = 
2.0 x 10 28 cm, eq. ([18]) becomes 



U = 2.5x10" 



10 ) 



10 40 ergs s- 1 sr- 1 



ergs cm 2 s 1 



for EC statistics. For synchrotron or SSC statistics 

Lf( £ 

10 ; VIO 42 ergs s" 1 sr" 1 



(19) 



f e 2.5x10" 



ergs cm 2 s , 



(20) 

so that a larger internal synchrotron power compared to 
Compton power is needed for the same vF v flux when 
viewing within the beam of the jet. 

In actual fitting of blazar statistical distributions, a 
degeneracy between T and l' e is found that may partially 
be removed by obtaining limits on the bul k Lorentz factor 
r from 77 attenuation argum ents (e.g. [Marasch i et al.l 
[1991 iDermer fe Gehrelslll995l) . 

For standard FSRQ blazar parameters, we take 



P 



3.4 , T = 10 , l' e = 10 40 ergs s" 1 



sr 



(21) 



and EC statistics. For standard BL blazar parameters, 
we take 



p = 3.0 , T = 4, l' e 
and syn/SSC statistics. 



10 42 ergs s~* sr" 1 , 



(22) 
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Fig. 3. — Model blazar formation rates (BFRs) used in this 
study. The curve labeled "cc" is the constant comoving rate, the 
curve "SFR," is the star fo rma tion rate from Hopkins & Beacom 
(2006), curve "IR,8" is eq. (12a I) with n = 8 that is used to fit the 
IR l um inosity density from IR luminous galaxies, curve "IR,4" is 
eq. (1 25 I t with n = 4 used to fit FSRQ data, and curve "BL" is eq. 
II26I I with a = 1.75 used to fit to the BL data. 

For two week (w 10 6 s) observations, « 50 flaring 
blazars were detected by EGRET above its vF v thresh- 
old flux fl hr = 2.4 x 10 -11 ergs cm -1 s _1 over w 1 year. 
Given the EGRET field-of-view (see Fig. 2), this implies 
a directional blazar flaring rate « 4 sr -1 Ms -1 , implying 
a fiducial rate density n com (z = 1) at redshift unity of 
about 



4x 10~ 6 s _1 sr" 



Hn 



4xlO bt W n com (z = 1) , (23) 



so that n com (z = 1) = 0.8 x 10~ events cm"^ s~ . 
Thus we expect that the blazar flaring rate density 
hcom = 10 -90 events cm~ 3 s _1 « 1 event Gpc~ 3 yr _1 . 
The calculations for the local rate density of FSRQs give 
values larger by one or two orders of magnitude. The 
discrepancy can be resolved by considering the increased 
number of sources necessary to offset the reduction in the 
detection rate by « T~ 2 for the beamed blazar emission, 
and the different rate density at z << 1 compared to 
z = 1 , which depends on the form of the BFR. 

3.5. Blazar Formation Rate Histories 

The BFR functions are used to describe the change 
in the rate density of blazars through cosmic time. This 
does not mean that the comoving density of supermassive 
black holes changes, but rather that their flaring rate 
changes or their assignment to the FSRQ or BL class 
changes with time as a consequnce of their definitions in 
terms of the optical emission line equivalent widths. 

We consider the following BFRs Tju{z), shown in Fig. 
[31 which give the blazar comoving rate density n com (z) — 
T,bi(z)n com in terms of the local comoving rate density 

1. Constant comoving rate density, so that S cc (z) = 
1. This form is employed for mathematical conve- 
nience. 

2. Comoving rate cx the blue/UV luminosity density, 
which is assumed to track the star formation rate 



G 



(SFR) of the universe. W e use the analytic form 
([Hopkins fc Beacomll200(l 



-'SFR 



(*) = 



1 + 6.78z 
1 + (z/3.3) 



5.2 



(24) 



Blazar activity could be related to the SFR if stellar 
activity provides fuel for the supermassive black 
hole engine, for example, from material driven off 
by starburst nurser ies encircling the nucleus (e.g., 
iHaiman et al.ll2004f ). 

3. Comoving rate oc sub-mm/far-IR lumino sity den- 
sity a ssociated with luminous IR galaxies ([Sanders! 
2004), which we fit using the analytic form 



^IR,n 



(*) = 



1 



(l + z)-" + 2- 



(25) 



We obtain a good fit to the data with n = 
8, as shown in Fig. 3. If IR-luminous galaxies 
are caused by galactic merger events, as is indi- 
cated by morphological a nd spectral evidence (e.g., 
iSanders fc Mirabell 19961 and references therein), 
this would connect blazars and the formation of 
supermassive black holes to galaxy collisions. Al- 
though related to supermassive black hole growth, 
the IR luminosity density does not, however, di- 
rectly measure the activity of supermassive black 
holes, because the IR radiation is a convolution of 
the photon luminosity which is then reprocessed 
through thick columns of material. Hence we have 
generalized the form with a single adjustable pa- 
rameter, n, that represents a range of BFR histo- 
ries. The forms of g and Y,m 4 are shown in 
Fig. 3. 

4. A BFR where the blazar flare rate density increases 
with cosmic time, which is found necessary to fit 
the BL data. The simple form considered is 



riBL(z) 



n B L{z = 1) 



(26) 



where a(> 0) is adjusted to fit the data. Because 
of the divergence in the rate density when z — ► 0, 
this BFR has to be normalized at z > 0. The total 
blazar flaring rate is, however, a well-defined value 
when a < 3. 

3.6. Diffuse Intensity 

The apparently diffuse intensity from the superposition 
of emissi ons from man y faint, unresolved blazars is given 
by (e.g.. lDermerll2006h 



2c 
A-kHq 



dz <b dQ! 



(1 



3 + Ov 

where e» = e z = (1 + z)e and a factor of 2 is again 
introduced for two-sided jet sources. The direction vector 
tt' defines the direction of the jet axis with respect to the 
observer direction. This expression applies to persistent 
blazar sources, with comoving emissivity 



d£~, 



dt*dV com dQ,' 



(28) 



When sources are detected above threshold flux f* hr , 
they no longer contribute to the EGRB, as they are iden- 
tified as a blazar source. This restricts eq. (|2T|) only to 
those blazars with VJ^t^" fd\ < fl hr , leading to the 
following expression for the diffuse radiation from unre- 
solved radio galaxies and black hole jet sources: 



(q - lJflo/JTfl Jo " ^/n m (i + z)3 + n A 

n com (z)l' e (z){[l - j 8min(l,/i)] 1 -« - (1 + Z?) 1 ^} . (29) 

The dependence of this expression on fl hr is carried by 
ft, given by eq. (fT5|) . 

Equating the EGRET two-week average fluxes with 
blazars that flare once in two weeks allows us to replace 
the comoving density of blazar AGN sources, given by 

/ \ / _ 1.2 x 10 6 s . . . _ 3 _! 
n com {z) (cm ) = — — x n com (z) cm s , 

L — p Z 

with the comoving rate density, n com (z), of blazar flares. 
This replacement should be accurate to a duty cycle fac- 
tor of order unity. Better studies based on GLAST obser- 
vations will reveal the flaring behavior of 7-ray blazars; 
note that the threshold flux fl hr ' G (t) for GLAST in the 
scanning mode is time-dependent when using eqs. (| 16|) 
and (TIT)) to fit the data. 

4. RESULTS 

Fig. 4 shows calculations of FSRQ blazar redshift dis- 
tributions using the standard FSRQ blazar parameter 
set, eq. (1211) . with differences from the standard parame- 
ters as labeled. In this calculation, the detection thresh- 
old is </>_8 = 15, corresponding to a vF v flux threshold 
sensitivity, from eq. |1}, °f 2.88 x 10~ 12 ergs cm -2 s _1 



a 4 10- 



BFR SFR (x 0.25) 




Standard Parameters: 
p=3.4, EC, BFR cc. r=10 
= 10 40 ergs s" 1 sr 1 , fli =15 

= 0.1,0.3,3,10 x10 4D ergs s 1 sr" 
(red dashed) 

r = 6,8,15,20 (blue dotted) 
= 10 39 ergs s" 1 sr" 1 



2 3 4 

Redshift z 



F = 6 / / 
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r= 6,8,15,20 (blue dotted) 




I' = 0.1,0.3,3,10 x10 4D ergs s" 1 sr" 1 




(red dashed) 
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Fig. 4. — Parameter study of the redshift distribution of FSRQs, 
showing the effects of different parameter choices on the differential 
redshift distributions (top) and cumulative redshift distributions 
(bottom). 
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at e — 200 for p — 3.4. Except for the light solid curves, 
which show results for a BFR given by the SFR history, 
eq. (pM)) , FSRQ blazars are assumed to emit blazar flares 
in random directions with a constant comoving rate den- 
sity, The local rate density, h com {z <C 1), is set equal to 
10" 90 cm" 3 s _1 . 

The effect of increasing r or l' e , of course, is to increase 
the distance from which blazars can be detected. Fig. 4a 
shows that the directional event rate per unit redshift is 
not greatly increased with increasing T factor, whereas 
this rate is increased with increasing l' e . The reason for 
this is that by increasing T, the emission is jetted into 
a smaller solid angle, rendering a smaller fraction of the 
blazars visible, though from a larger distance. Neverthe- 
less, the cumulative redshift distributions for different 
combinations of T and l' e can be very similar, as seen in 
Fig. 4b. When r > 1, these two quantities enter into 
the threshold cosine angle p, according to the combina- 
tion l f J~' q /T (eq. [j"5j). reflecting the degeneracy of the 
results as a function of these two quantities. 

By changing the BFR function, the distribution of 
blazars with redshift can be adjusted to better agree with 
the observations. As is made clear by Fig. 4b, the effect 
of using the SFR function rather than the constant co- 
moving rate is to decrease the number of low-redshift 
blazars, and to have more blazars detected at 1 < z < 3. 

Even though it is relatively simple to find a parameter 
set and BFR that gives a good fit to the FSRQ redshift 
distribution, it is considerably less simple to find a pa- 
rameter set that gives a good fit to the joint redshift and 
size distribution. One difficulty in fitting the size distri- 
bution originates from sample incompleteness near the 
EGRET threshold. The value 4>s = 15 applies to two- 
week observations of on- axis EGRET sources (Appendix 
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0.2 - 
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/ FSRQ 




f r=10, p = 3.4 




f I' = 2.5x1 39 ergs s' 1 sr' 1 




A? EC, BFR IR,4 




T » =25 
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Redshift z 
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I' = 2.5x10 39 ergss' 1 sr' 1 

EC, BFR IR,4 
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Peak flux if [1 0" a (ph> 1 00 MeV) cm" z s" 1 ] 

Fig. 5. — Smooth curves give model fits to the cumulative red- 
shift and size distributions of FSRQ EGRET data, shown by the 
histograms in the top and bottom panels, respectively. Parameters 
for the fit are shown in the legend to the figure. 



A). Reduction in the effective area for off-axis sources 
means that the effective threshold where the EGRET 
sample is complete is larger than <f>-s = 15. As a con- 
sequence, in the fitting of the data for FSRQs and BLs, 
the flux threshold was changed to <j>^s = 25 where the 
EGRET sample evidently no longer suffers from incom- 
pleteness. After adjusting the parameters and BFR his- 
tories, an acceptable fit to the FSRQ was obtained, as 
shown in Fig. 5. 

This fit uses BFR IR,4 shown in Fig. 3, with param- 
eter values shown in the figure caption. The value of 
p = 3.4 was chosen to to agre e with the mean FSRQ 
photon spectral index a p h = 2.2 (Mukher iee et al.lll997| ). 
The values of T and l' e are not unique due to the de- 
generacy mentioned above, but the choice of T — 10 is 
suggested by superlumi nal radio observations of FSRQs 
(|Urrv fc Padovanilll995h . As can be seen, the model fits 
the cumulative size distribution at <^)_8 > 25 but not at 
lower values of 4>s due to sample incompleteness. The 
FSRQ model gives a statistically acceptable fit to the 
distributions, noting that the Kolmogorov-Smirnov one- 
sided statistic for 46 sources (42 sources with = 25.) 
is 0.16 and 0.22 at the 90% and and 99% confidence 
level, respectivly. The model fit, does, however, show 
a slight deficit of blazars at 1.5 < z < 2 compared to the 
data. Moreover, the model somewhat overproduces the 
number of very bright FSRQs. Both the high-redshift 
deficit and overproduction of the brightest FSRQs could 
be alleviated by tuning the BFR to increase even faster 
than given by the model BFR, inasmuch as the bright- 
est sources are generally found at lower redshifts. Given 
the statistically acceptable fit, further fine-tuning would, 
however, introduce additional parameters that are not 
well constrained. 

Even after searching over a wide range of BFRs for 
the BLs, it was not found possible to obtain acceptable 
fits to the EGRET data for the joint redshift and size 
distributions of BLs. The difficulty was the very narrow 
range of peak fluxes, spanning less than a factor sa 5 
from the dimmest to the brightest values (compared to a 
factor ps 20 for the FSRQs), and the requirement to have 
the same threshold integral photon flux </>_g = 25. The 
crux of the problem is that the blazar size distribution is 
steeper than —3/2 (Fig. 2). A resolution of this problem 
was to introduce luminosity evolution of the BLs such 
that they dimmed with increasing time. At the same 
time, the comoving rate density of BLs increases with 
time, allowing for a large number of nearby BLs to be 
observed so that the redshift distribution of BLs could 
be fit. The effect of dimming luminosity and increasing 
rate density with time makes blazar fluxes over a narrow 
range of values. 

Fig. 6 shows the result of this procedure, using eq. 
i{26"|) to describe the BFR of BLs, with a = 1.75, and 
luminosity evolution described by l' e oc z 195 . We choose 
p = 2 to agree with the mean BL ph oton spectral in- 
dex a p h = 2.0 (jMukherjee et al.l 119971 ) . which implies a 
threshold vF v EGRET flux of ff r = 4 x 10" 11 ergs 
cm~ 2 s _1 for (j> = s — 25. Other parameters of this model 
are given in the caption to Fig. 6. The model gives a 
statistically acceptable fit, noting that in this case, the 
Kolmogorov-Smirnov one-sided statistic for 14 sources 
(12 sources for 0_ 8 = 25.) is 0.275 and 0.39 at the 90% 



and and 99% confidence level, respectivly. 

The sample of FSRQs and BLs with (fi-g > 25 suffers 
much less from sample incompleteness near threshold. 
This flux-limited sample has 4 fewer FSRQ and 2 fewer 
BL sources, and is therefore somewhat less constraining 
to the model fits. Although the fits to the size distribu- 
tion using this smaller sample would be improved near 
threshold, the modified redshift distribution is not sig- 
nificantly changed, and the use of this sample does not 
change the conclusions of this study. 

5. PREDICTIONS 

Fig. Q] shows the fits to the FSRQ and BL redshift dis- 
tributions implied by the models discussed in the previ- 
ous section. The factors used to normalize to the flaring 
rates give a local FSRQ flare rate density equal to 

h FS RQ = 5.66 x 1CT 88 cm' 3 s" 1 = 
17 Gpc" 3 Ms" 1 = 1.7 x 1(T 8 Mpc~ 3 Ms -1 (30) 

for FSRQs, and a comoving BL flare rate density at z = 1 
equal to 

u B l{z = 1) = 1.23 x 1(T 88 cm -3 s" 1 = 

3.6 Gpc" 3 Ms" 1 S 3.6 x 1CT 9 Mpc~ 3 Ms" 1 (31) 

for BLs. The predicted size distributions of FSRQs and 
BLs extrapoloated to small values of 0_8 implied by the 
model are shown by the solid curves in Fig. O Also 
shown is the one-year sensitivity of GLAST in the scan- 
ning mode, which is at the level of « 0.4 x 10~ 8 ph(> 100 
MeV) cm" 2 s -1 . The model fit predicts that GLAST will 
detect « 700 FSRQs and FR2 radio galaxies, and w 160 
BLs and FR1 radio galaxies after one year of observa- 
tion, and an additional ss 200 dim FSRQs/FR2 galaxies 




Peak flux t, [10" B (ph> 100 MeV) cm" z s" 1 ] 

Fig. 6. — Smooth curves give model fits to the cumulative redshift 
and size distributions of BL EGRET data, shown by the histograms 
in the top and bottom panels, respectively. Parameters for the fit 
are shown in the legend to the figure. The rate-density dependence 
of the BL Lac flares is oc 2~ 1,75 . 



and w 50 more BLs/FRl galaxies after 3 more years of 
observation. These numbers should be taken with an 
estimated uncertainty less than a factor-of-2, consider- 
ing our lack of knowledge of blazar activity at z > 3. 
The predicted numbers would represent a lower limit if 
there are blazar subclasses yet to be discovered that are 
numerous. Broadening of the luminosity and T distribu- 
tions might increase the total number, though by only 
factor of order unity. Additional classes of low luminos- 
ity blazars that were not detected with EGRET could 
be discovered with GLAST, and these would produce an 
upturn in the size distribution at low fluxes. The number 
of such sources are, of course, limited by the level of the 
EGRB. 

Fig. 7 shows the predicted FSRQ and BL blazar red- 
shift distributions for GLAST. The distribution pre- 
dicted for EGRET is shown by the dashed curves, corre- 
sponding to the result shown in Fig. 1, and the predicted 
distributions for GLAST sensitivities of 1,3,10,30, and 
100 x the EGRET sensitivity of 4>s = 15 are shown by 
the solid curves. 

Fig. 8 shows the contribution to the E GRB as reported 
by th e EGRET team (filled data points: ISreekumar et alJ 
1998), and as derived from the EGRET data using 
the GA LPROP cosmic-ray propa gation code (open data 
points; IStrong et all |2000L 12004( 1 . The boxes give the 
contributions to the EGRB from unresolved FSRQs and 
BLs with flux values between <^_g = 25, below which 
a two-week observation is incomplete, and </>_g = 12.5. 
This latter value is a factor of two less than the brighter 
value, and is chosen because EGRET had an effective 
lifetime of « 4 years in terms of its capability during the 
first year, due to limited amount of spark chamber gas. 
This would produce a factor w 2 better sensitivity before 
sample incompleteness dominates. Nonuniformity in ex- 
posure and angle to the normal of EGRET suggests that 
the threshold is in fact quite broad. 

The steeper FSRQ and flatter BL spectra in Fig. 8 
arise from using p = 3.4 for FSRQs andp = 3 for BL Lacs 



300 




Redshift z 

Fig. 7. — Histograms show the redshift distributions of FSRQs 
and BLs discovered with EGRET. Dotted curves show the model 
fits to the distributions, using detection characteristics typical of 
EGRET (and 4>-g = 25; see Fig. 1). Solid curves show the pre- 
dicted FSRQ and BL redshift distributions using the parameters 
used to model the EGRET data, though with detection character- 
istics of GLAST. For the GLAST predictions, the observing photon 
energy is set equal to 1 GeV (e f*s 2000), and the detection sensi- 
tivity is 1, 3, 10, 30, and 100X better than EGRET, using the 
4>S = 15 as the flux threshold for EGRET. 
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Fig. 8. — The intensity of the isotropic EG RB measured wit h 
EGRET is shown by the filled data points (Sreckumar ct al,l ll998T) , 
and the EGRB after subtraction of a component of q uasi-isotropic 
galac tic 7 radiation is shown by the open data points (jStrong et al.l 
|2004f ). The dashed box gives the diffuse intensity from unresolved 
FSRQs, and the dotted box give the diffuse intensity from unre- 
solved BLs, and the curve labeled "Total AGN" is the sum of the 
two. The boundary of the box represents the flux limit below which 
blazars arc assumed to be unresolved, and correspond to the range 
for 4>—% = 25 and </>_g = 12.5. 

(note also the weak curvature in the calculated FSRQ 
background spectrum). The model results imply that 
the total AGN (BL+FSRQ) con t ributi on to the EGRB 
as measured by ISreekumar et alj (|1998l ) is at the level of 
w 12 - 20%, and has a spectral index remarkably simi- 
lar to the index measured in the EGRET analysis. The 
FSRQ contribution to the EGRET 7-ray background at 
1 GeV is estimated at the level of ss 10 - 15%, and the 
BL contribution is sa 2% - 4%. This co ntribution also 
includ es FR1 radio galaxy, estimated by IStawarz et al.l 
(2006) as contributing at the « 1% level (« 6 x 10" 3 
keV cm~ 2 s _1 sr -1 . The other sources classes that con- 
tribute to the residual EGRB are considered elsewhere 
(jDermerl [20061 . Note that no hard component BL Lac 
contribution to the EGRB is calculated here, though one 
must exist given meas urements of hard-G e V spectrum 
BL L ac objects (e.g., iKataoka et al.1 [T999t iPetrv et"afl 
2000, for Mrk 501, which is not in our sample). Spec- 
tral hardening during flaring is also not considered here. 
Even so, this calculation implies that there is room for 
a number of other sources and/or source classes to make 
up the EGRB. 

6. DISCUSSION 

We have employed a simple physical model to fit pop- 
ulation statistics of EGRET 7-ray blazar data and make 
predictions for the data anticipated from GLAST. By 
characterizing all FSRQs and BLs in terms of single val- 
ues of r and l' e , the model contains the fewest num- 
ber of free parameters while still agreeing with the ba- 
sic physical understanding of blazars as a two-sided jet 
of collimated relativistic plasma outflow with the radi- 
ation Doppler boosted by the motion of the jet. Other 
than blazar rate density evolution, described in terms 



of the BFR function, only the parameters T, the proper 
directional power V el the power law index p and beam- 
ing statistic q define a model. The values of p and q 
are known with some certainty from observations. The 
model is further simplified in that the cumulative dis- 
tributions are degenerate in the quantity V^ q jY when 

r > 1. 

We have used r = 10 for FSRQs and T = 4 
for BLs. These values of T are generally consis- 
tent with Lorentz factors inferred from fits to the lo- 
cal luminosity functions of radio galaxies and blazars 
(|Urrv fc Padovanil I1995D and measurements of super- 
luminal motions (|Vermeulen fc Cohenl 1 19941 ) . The su- 
perluminal motion observations indicate that BLs have 
smaller Doppler factors on average than FSRQs, though 
in both cases the spr ead in Doppler factors is rather large 
(|Jorstad et alJ2001f ). Good fits to the FSRQ data are ob- 
tained with parameters T = 10, l' e = 2.5 x 10 39 ergs s _1 
sr , p = 3.4, and the EC beaming statistic. A BFR 
function that increases by a factor 10 - 20 from z <C 1 to 
z « 1 - 2 was required to jointly fit the FSRQ redshift 
and size distribution, but a good fit was obtained even 
without luminosity evolution. 

In contrast, it was not found possible to obtain a good 
fit to the BL redshift and size distributions without also 
invoking luminosity evolution. The BL size distribution 
declines more rapidly than a uniform Euclidean size dis- 



tribution 2V(> 



oc 



-3/2 



as can be seen in Fig. 2. 



Any luminosity function, whether a kinematic luminos- 
ity function produced by the randomly oriented blazar 
jets, or a luminosity function given in terms of a range of 
jet powers, only flattens the size distribution. A feasible 
solution to this dilemma was to have the BL luminosity 
dim with time while the rate density of BL flares in- 
creases with time. A model that provides a reasonable 
fit totheBL data has T = 4, l' e (z) = 6xl0 4 V' 9S ergss" 1 
sr -1 , a rate density evolution of BLs cx z~ 9 / 4 , p = 3.0, 
and the synchrotron/SSC beaming statistic. 

Before proceeding to the implications of this analysis, 
a few points should be made. 

1. The tightly constrained parameter space in this 
simplified blazar model leads to fairly robust pre- 
dictions for the number and intensity of blazars, 
and we think that the number of blazars and radio 
galaxies cannot be more than a factor-of-two in er- 
ror, and then in the direction of more sources and 
greater AGN background 7 rays. A more precise 
error analysis must be considered in further study 
of this model. 

2. Luminosity evolution was introduced in the BL fit- 
ting, but not the FSRQ fitting. The blazar's power 
evolves with cosmic times, reflecting the growth 
of the supermassive black hole (including binary 
black-hole merger episodes) and its intermittent 
fueling. The jet power is related to the accre- 
tion power in ways not fully understood, and is 
likely also dependent on black-hole spin. Modeling 
FSRQ growth and fueling with only a BFR function 
is possible with the EGRET data, but may well 
fail when modeling the larger and more meaningful 
GLAST data set. The next step in the analysis ef- 
fort is to allow luminosity evolution of both FSRQs 
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and BLs based on a physical model for supermas- 
sive black-hole growth and fueling, including duty 
factors on small and large timescales. 

3. The quality of the fits, tho ugh acceptable , is no t 
as good as in the study of iMiicke fe Pohll ([2000). 
This indicates that a range of luminosities or in- 
jection energies may be required to obtain the best 
fits for the FSRQs, and that the form of the 7- 
ray luminosity function may be well represented 
by the luminosity function of FR1 and FR2 radio 
galaxies. Note also th at the BL sample used in the 
IMiicke fc Pohll (1200(1 appears smaller and less con- 
straining than the BL sample used here. The size 
distributions of BLs measured with GLAST will 
show whether the BL sample used in this study is 
representative. 

6.1. Blazar evolution 

The solution presented here that fits both the steep 
BL size distribution and the BL redshift distribution was 
to have joint positive luminosity evolution and negative 
source density evolution. From a simple size-distribution 
analysis in the Euclidean limit, the flux ((/>) dependence 
of the integral size distribution goes as 

JV(> 4>) CX 0-(3-a„)/(2-a,) ^ 

where l' e oc z ai and n(z) oc z _a ™, and a n = for a 
uniform distribution of sources. A flatter (steeper) size 
distribution occurs for negative (positive) luminosity evo- 
lution with ai < (>)0. Arranging to fit not only the BL 
size distribution with slope steeper than the uniform Eu- 
clidean value of —3/2 requires positive luminosity evolu- 
tion, and also to fit the measured BL redshift distribution 
required negative source rate-density evolution, leading 
after trial and error to the result a; = 1.95 and a n — 9/4. 
Recalling the definition of (V/V max ) as the normalized 

sum over the test statistic (4>i/<j>thr)~ 3 ^ 2 , one sees that 
size distributions steeper than —3/2 have a larger num- 
ber of faint detections near threshold than for a uniform 
Euclidean distributions and thus have (V/V max ) > 0.5; 
size distributions flatter than —3/2 have a deficit of faint 
detections near threshold and thus have (V/V max ) < 0.5. 

This EGRET BL analysis indicates that (V/V max ) > 
0.5, similar to the case for bright 1 Jy radio sources or 
the brightest BL Lac objects selected at X-ray energies 
(jGiommi et al.l l200ll ). From Fig. 2, we predict that 
{V/Vmax) for BLs will shift to < 0.5 within a few months 
of the start of GLAST measurements, as fainter BL Lac 
objects and FR1 radio galaxies will be lost due to cosmo- 
logical effects and the small number of high-redshift BLs. 
Our analysis leads, however, to very different conclusions 
than obtained previously when considering only the size 
distribution and (V/V max ) value from various BL surveys 
at radio and X-ray energies. For example, a flatter than 
-3/2 size distribution is found in the ROSAT (0.5 - 2.0 
keV) all-sky survey (RASS) of over 30 BL Lac objects 
complete to ft (0.5 - 2.0 keV) > 8 x 10~ 13 ergs cm~ 2 s" 1 
(|Bade et al.|[l998l ). corresponding to a deficit of sources 
near threshold, or negative B L luminosity evolution . The 
same conclusion is found by iGiommi et alJ ([1999) , who 
assemble a sample of 155 BL Lacs by cross correlating 
the RASS survey results with radio sources from the 



NRAO VLA Sky Survey. The resulting radio luminos- 
ity size distribution of these X-ray selecte d BL Lac ob- 
j ects a lso exhibit a negative BL evolution. iRector et al.l 
(|2000f ) measured redshifts and analyzed sources in the 
Einstein 0.3 - 3.5 Medium Sensitivity Survey, and ob- 
tained (V/Vmax) — 0.4, again showing a deficit of near 
threshold sources compared to the uniform, Euclidean 
expectation. 

It is meaningful to interpret a flattening of the size 
distribution compared to a —3/ 2 slope as a reduc- 
tion in luminosity at early times (|Padovani fc Giomrtril 
1995T). Here howeve r we additionally require, as do 
Miicke fe Pohl (|2000f ). that the redshift distribution is 
also accurately fit with the same model that fits the size 
distribution. Unlike the lGiommi et all (|2006f ) calculation 
of the EGRB from 7-ray blazars, no radio/7-ray connec- 
tion is assumed throughout this analysis (though radio 
data is used to guide the choices of T). The results from 
joint fitting of EGRET data are severely constraining: 
the comoving rate density of FSRQ declines by a fac- 
tor > 10 at z <C 1 compared to z < 1 - 2, while the 
rate density of BLs rapidly increases for z <C 1 and, at 
the same time, the mean luminosity of BL flares declines 
with time. 

This behavior is in striki n g accord with the s c enario 
(jBottcher fc Dermerl l2002t ICavaliere fc D'Elial l2002f ) 
linking BLs to FSRQs through an evolution in declin- 
ing accretion rates and incre asing black hole ma s ses. I n 
the particular formulation of lBottcher fc Dermerl ([2002) , 
the reduction in the amount of gas and dust that both 
fuels the central black hole engine and scatters accretion- 
disk radiation causes a transformation of the FSRQs 
into the BLs in terms of spectral properties, thereby ex- 
plaining the blazar sequence correlating vF v peak syn- 
chrotron frequencies and apparent synch rotron power 
(jSambruna et al]ll996t iFossati et aLlll998l ). and puts on 
a more physical basis the correlations between the peak 
synchrotron frequency, peak electron Lorentz factor, and 
inj ection and external ra diation compactnesses proposed 
bylGh isellini et al. (1991) to explain the blazar sequence. 
ICavaliere fc D'Elial (|2002h argue that the mean luminos- 
ity of the BLs does not change much over several gi- 
gayears, and that the flattened size distribution of BLs 
at low fluxes stem from negative density evolution, with 
spectral differences associated with transition from ac- 
cretion power to a component from the black-hole spin. 

A clearcut prediction of this scenario is that the mean 
masses of black holes in BLs is larger than that in FSRQs 
at the same epoch. The minimum variability timescale 
of blazar flares should be proportional to black hole mass 
(§A.5), but the shorter timescale flaring behavior of BL 
Lacs measured with Whipple and HESS compared with 
the minimum variability timescale of FSRQs measured 
with EGRET is opposite to the expected behavior, but 
may only reflect sensitivity limitations of EGRET. Prob- 
ably there is flaring on every timescale, and a power spec- 
tral density analysis is required to see if there is a size 
scale where temporal power declines, corresponding to 
the black hole mass. This question will obviously be 
s ubject to extensive investig ation with GLAST. 

iGhisellini fc Celottil (|200lh infer the black hole masses 
of FR1 and FR2 radio galaxies on the basis of an ex- 
pression for the central black hole mass in terms of host 
galaxy absolute R-band magnitudes. Relating jet power 
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to radio luminosity suggests that FR1 galaxies are low 
Eddington-ratio (< 1CP 2 ) and FR2 galaxies are mod- 
erate Eddington-ratio (0.01 - 0.1) sources. There is a 
tendency for the FR1 galaxies to have larger black-hole 
masses than FR2 galaxies, but this effect is too weak to 
provide an unambiguous test of the FSRQ — > BL evo- 
lutionary scenario. 

It is unphysical to extend the form of the BL rate den- 
sity evolution, oc z~ 9 / 4 , to arbitrarily small redshifts. 
The comoving density of FSRQs forming at high red- 
shifts should equal the comoving density of the BLs at 
low redshifts in the picture of FSRQ — > BL evolution 
considered here. The BL density will saturate at zu 
when all FSRQs have converted into BLs. If the duty 
cycle of the two sources classes is similar, this occurs 
at tibl{zu) — riFSRQ{ z ^ !)■ Using the derived values 
in eqs. ([30]) and (f3"Tjl and the BFR functions, this takes 
place at Zu = 0.1. Modifying the BL BFR to be constant 
at z < zu has minimal effects on the results, but should 
be considered in the next step of the analysis. 

6.2. Space density of blazars 

The local density of blazar sources given from this 
analysis can be compared with the space density of 
blazar host galaxies, namely early-type elliptical galax- 
ies. Fro m a K-band survey of bright galaxies with 
z < 0.4, iHuang et al.l (|2003[ ) calculate a local galaxy 
density of 0.0048 ± 0.001 Mpc -3 for galaxies of all types 
with h = 0.72. The density of early-type galaxies is 
« 10% of the total, or w 5 x 10~ 4 Mpc~ 3 . In compari- 
son, the comoving densities of the sources of FSRQ and 
BL flares are Ufsrq(z « 1) ~ 2 x 10~ 8 Mpc~ 3 and 
nBh{z = 1) — 2 x 10~ 9 Mpc -3 , assuming a duty factor 
equal to unity (a duty factor less than unity implies a 
proportionately smaller source density). 

The comoving density of sources at the present epoch 
is, in the scenario where FSRQs evolve into BLs, equal to 
~ Z FSRQ (z > l)n FSRQ (z < 1) = l7n FSRQ (z < 1) = 
4 x 10 -7 Mpc -3 . This is about three orders of magnitude 
smaller than the elliptical galaxy density. Either only a 
small fraction of ellipticals host radio galaxies, or the 
duty cycle of elliptical galaxies that harbor radio and 7- 
ray jets is a small fraction of the total lifetime of the 
source (see lHaiman et alll2004 for more speculations). 

7. SUMMARY AND CONCLUSIONS 

A method to analyze the population statistics of 7- 
ray blazars solely f rom the 7-ray data w as considered in 
this paper (see also lMucke fc Pohll2000l ). By performing 
the 7-ray population study independent of other wave- 
bands, the reliability of this method to other methods 
that invoke a radio/7-ray correlation can be tested, and 
the physical reasons for differences in radio galaxy and 
blazar populations selected from observations at different 
wavebands can be explored. 

Although a radio/7-ray correlation was avoided, guid- 
ance to assign values of T was taken from radio obser- 
vations. Measurements of apparent super luminal motion 
from radio observations on the scale of as 0.1 - 1 pc from 
the black hole imply T factors. The 7 rays may origi- 
nate from within hundreds to thousands of Schwarzschild 
radii of the supermassive black hole, so that the T values 
derived from radio observations may not be appropri- 
ate to the analysis of 7-ray data. Radio observations 



suggest that 7-ray flares could originate f rom the same 
physic al scale as the superluminal blobs (|Jorstad et al.l 
l2001al) . though spectral modeling of some X-ray se- 
lected BLs such as Mrk 4 21 and Mrk 501 suggest l arger 
Dopp ler factors ~ 50 ([K rawczyns ki et al.l l200ll but 
see iGeorganopoulos fc Kazanasl (|2003l ): iGhisellini et al.l 
(|2005f ) for possible resolutions to this puzzle) than in- 
ferred from superluminal motion observations. Lower 
limits to the values of T for BLs and FSRQs will be 
inferred from 77 attenuation arguments applied to the 
GLAST observations. This will not only provide better 
values to use for modeling the population statistics of 
7-ray blazars, but will help break the parameter degen- 
eracy in the modeling. 

By avoiding an underlying radio/7-ray assumption, 
only a very simplified blazar model could be investi- 
gated. This model nevertheless contains the essential 
blazar physics, and can be simply generalized to include 
a range of T factors and an evolving, broadened lumi- 
nosity function. For the analysis of the EGRET data, 
introducing distributions in T or l' e would allow too large 
a parameter space to explore without introducing some 
radio/7-ray connections, so we considered fixed values of 
T = 10 for FSRQs and T = 4 for BLs, and found val- 
ues of luminosity l' e that gave reasonable fits to the data. 
Because mean spectral indices were assigned from 7-ray 
observations, and the beaming factor statistic was taken 
from blazar physics, only the value of l' e , the BFR func- 
tion and, when necessary, the redshift dependence of l' e , 
were varied in order to fit the EGRET redshift and size 
distributions. 

Within these constraints, the EGRET data for FSRQs 
was fit with a BFR function that had ~ 10 x more sources 
at z = 1 than at present, which could be related ei- 
ther to a star-formation rate function or an evolutionary 
behavior proportional to the far IR/sub -millimeter lu- 
minosity density related to IR radiation (jSanderj 120041 : 
iBlazeiowski et al.l I2000D . The EGRET data for BLs 
could not be fit by only modifying the form of the BFR, 
and luminosity evolution was also required. The result- 
ing solution — that the density of BLs increase and their 
mean jet powers dim with time — is in accord with the 
scenario where FSRQs evolve from BLs. 

The analysis implies that the contribution of unre- 
solved blazars to the e xtragalactic 7-ray back ground 
measured with EGRET (jSreekumar et alj|1998T) at the 
level of « 20% at 1 GeV, leaving room for various source 
classes or additional types of blazars, such as hard GeV- 
spectrum blazars, to which EGRET was not very sensi- 
tive. Approximately 1000 blazars are predicted from this 
analysis to be discovered with GLAST during its first 
year of operation, but this could be an underestimate for 
the reason just stated, though not by more than a fac- 
tor of 2. As also observed with EGRET, the BL/FR1 
7-ray sources are predicted to have a much smaller mean 
redshift than the FSRQ/FR2 sources (Fig. 7). The pre- 
dictions for the high redshift (z > 3 - 5) blazar popula- 
tion is very uncertain because of the very few high-z 7- 
ray blazars observed with EGRET. If blazars only make 
~ 20% of the EGRB at 1 GeV, then star-forming galax- 
ies, starburst galaxies, and cluster of galaxies are likely 
to contribute a significant fraction of the difference, so we 
can expect GLAST to discover new 7-ray source classes. 

With the much larger data set from GLAST, more 
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detailed analyses will allow various effects related to 
blazar spectra and flaring and the evolutionary connec- 
tions between various classes of cosmological black-hole 
jet sources to be examined in much greater detail. 
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APPENDIX 

SENSITIVITY OF EGRET AND GLAST TO BLAZAR FLARES 



We follow the approach of iDermer &: Dingu s (2003). The properties of a 7-ray imaging spark chamber or silicon 
tracker detector depends on its shower pattern, assumed to be described by a Gaussian with 68% containment angle 
9 and energy-dependent angular point spread function (psf) described by 



fpsf 



9u 



where the photon energy in units of 100 MeV is 



E 



£1 



and E 



100 



100 MeV 



For EGRET, 9 ^ 5.7°/57.3° ^ 0.1 and w = 1/2. For GLAST, 9 3.5°/57.3° 0.06 and w = 2/3. 

The effective area A(E, 9, <fi) depends on photon energy E, angle 9 from zenith, and azimuthal angle 
with respect to the zenith angle, according to 



A(E, 9, < 



A{u) £ A u% 



(Al) 

(A2) 

measured 
(A3) 



where A is the effective area averaged over the field of view, defined as the opening solid angle within which the 
effective area is 50% of the on-axis FOV. The FOV of EGRET is = 4tt/24 S 0.5 sr, and the FOV of GLAST is 
= 47r/6 = 2 sr. The parameters describing the effective area Aq at 100 MeV and index ao are, respectively, = 1200 
cm 2 and a = for EGRET, and S 6200 cm 2 and a a = 0.16 for GLAST. 3 
The significance to detect a signal at the no~ level is given by 



S 



y/aS + (1 + a)B 



(A4) 



(|Li fc Ma 1983, eq. 9), where S is the number of source counts, B is the number of background counts, and a — t on /t f / 
is the ratio of on-source to off-source observing times. If the background in precisely known, a — ► 0, and 



S_ 
B 



(A5) 

In the limit a <C 1 and n 3> 1, eqs. (1) and (2) overestimat e the fluctuatio n probability compared to a maximum 
likelihood expression confirmed by Monte Carlo simulations f|Li Sz Malll983f ). so eq. (|A5[) should be a conservative 
expression for the detection significance. 

Background Counts 

We assume that the source location is precisely known. The number of background photons within solid element 
AQ(E) centered in the direction f2, and with energies > E\ observed during an observing time At is 



At 



B(>Ex)^ dt dE ACl(E) A[E,9(t),cf)(t)}<f>B{E,n) 



(A6) 



Ei 



The background photon flux per steradian, &b(E, f2) = dN 7 / dAdtdfldE , is assumed to be time-independent; thus this 
treatment excludes passage through the South Atlantic Anomaly. This expression also applies to time-independent 
diffuse backgrounds, not time-variable point sources or variable backgrounds. The effective area A[E, 9(t), <j>(t)] of the 
telescope at energy E changes for a source at time- varying angles 9 and <f> with respect to the telescope z-axis, which 
change with time due to rocking or slewing or Earth occultation. These effects are taken into account with an exposure 
factor X. Note that eq. (|A6[) assumes that each photon count can be precisely assigned an energy E and direction 



3 These satisfy the GLAST Science Requirements Document; the actual performance is given at the GLAST websites. Note that the 
estimate requires an average over the FOV. 
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The apparently diffuse isotropic 7-ray background spectrum is independent of f2 and is given by ISreekumar et al.l 
(1998) as 

* X (E) = k x u- aB , (A7) 

where k x — 1.73 ± 0.08 x 10~ 7 ph (cm 2 -s-sr-MeV) _1 and as = 2.10 ± 0.03. This expression is a valid description of 
the minimum 7-ray background i rregardless of whet her the it is extragalactic, or contains a quasi-isotropic galactic or 
heliospheric emission component (|Strong et al.ll2004l ). Thus the number of background counts with energy E observed 
over the real time At — t w k weeks is 

= XAtE 100 n9 2 A k x 1+ao _ 2 _ aB 
v ; 2w + a B -l-a ' 

Source Counts 

For a point source, the number of source counts with energy > E is given by 

r At 

dE A[E,e(t),(/)(t)] <j>s{E,t) , (A9) 



plAZ POO 

S(> E,) = / 7 / dt 

J0 J Ei 



where <p s (E,t) is the source photon flux (ph cm 2 s 1 E 1 ) and we assume that the point spread function is defined 
as the f 1 = 68% containment radius (see lThompsonl [l986) . Thus 

s 10-%X s At(a ph - 1)^_ 8 ul+ao _ aph 
otph — ao — 1 

and a p ft, is the photon number index (commonly denoted T). The factor X s accounts for the on-source exposure, and 
is generally equal to the background exposure factor X, though it may differ if one allows X s to account for the source 
duty cycle, of if X accounts for variable background. The quantity normalizes the source flux 

ME) = 10 ' 8( ^ = 1)0 - 8 «-•» (All) 

to 10~ 8 ph(E > 100 Me V) cm" 2 s" 1 , which is the unit quoted by the EGRET team in the Third EGRET catalog 
(jHartman et al.llT999l ) for the two-week average fluxes. 

For EGRET, S(> u) = 4.9u~ 11 0_8^«)fe > 5ri5 S , where a source detection is assumed to require at least 5 counts 
(ri5 S = 1). From the signal limit, a minimum time 

u <*ph-l 2u 1A , 

% —^T n 5 S -— , (A12) 

is needed for EGRET to detect sources in the EGRET FOV as faint as 4>-$ (assuming that the blazar remains at 
the level of <j)_g during the entire viewing period). Earth occultation means X s f» 1/2. The final expression in eq. 
(I A12|) applies when a p h = 2.1, the index of the isotropic 7-ray background. Here we suppose that blazars make a 
large fraction of the EGRB, and so must have a significant fraction of bright blazar sources emitting with a p h = 2.3 
between w 100 MeV and 10 GeV. 
For GLAST in the scanning mode, the condition for the detection of 5 counts is 

*Sfc £ — — r y a, " 5s - °- 85 -7— ' ( A13 ) 

V a ph - 1 / A!/ 5 0_ 8 0_ 8 

where we have adopted a source occultation factor X s = O.2X1/5. The two estimates are similar at 100 MeV, because 
the factor 6-8 effective area advantage of GLAST over EGRET is reduced by GLAST in the scanning mode compared 
with the few occasions when the source is a pointing target for EGRET. 

Signal to Background 
From eqs. (|A5[) . (|A8|1 and (|A10[) . the significance to detect a blazar at the na level is 



where 



and 



S(> u) lO- 8 / 7 A,V^Al0„ 8 

n(> u) = ; = — (p c u Vx , (A14) 

VB(> u) VXnk x E W0 6 

= ^ ~ 1} : V2^ + a B - 1 - a , (A15) 
a P h —ao — 1 

l + a + a B 

<Px = w — oiph H • (A16) 
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The energy dependence of the background-limited detections of blazars is defined by the index 

(2.05-a ph , EGRET 
Vx = { ' , (A17) 

[ 2.30 - a ph , GLAST 

and it is of interest that EGRET and GLAST are more sensitive to hard-spectrum blazars (a p h < 2 - 2.3) at higher 
photon energies. This sensitivity improves until either there is a spectral break or softening, or the signal runs out. 
The background limit on EGRET, with X = X s , is 

n E (> u) = 0.36V^W-8U 2 '° 5 ~ Qph > 5n 5(T . (A18) 
Taking X = 0.5, due to Earth occultation, the number of weeks of direct pointing required for EGRET to detect a 
source with flux at the level <^>_g at the 5n 5rT a significance level is 

& £ ? " 2(Qph " 2 -° 5) "L ■ (A19) 

Comparison with eq. (IA12|) shows, as is well known, that signal detection with EGRET is background rather than 
signal limited. For u = 1, as used in the analysis reported in the EGRET catalogs, the limiting sensitivity to detect a 
strong source in a two-week pointing is 

0-8 > 14 • 

This estimate agrees well with detection data, as seen from Fig. 2. 
For GLAST, the background limit translates into the condition 

3.6n 

»-8 - 



> ^!*o_ u °p*-2-3 t ( A 20) 



yXtwk 

so that the number of weeks of scanning required for GLAST to detect a source with flux level 0„8 at the ^n^a 
significance level is 

& £ y ,2 5 - "3" «" 2/5 ■ ( A21 ) 

Thus, in around a year in the scanning mode, GLAST will detect sources at the level 0_8 ~ 1 when integrating above 
100 MeV (u — 1). This is not however the best detection strategy. 

Optimal Source Detection with GLAST 

A better analysis considers optimal energy for source detection from both signal and noise limits. Equating eqs. (|A13[) 
and (|A21|) for the minimum time needed to detect sources at the 5a, 5 count limit assuming p = 3.2 or a v = —0.1, 
gives the best lower photon energy above which to integrate, namely 

u = 20 4>Zl /A , or E = 2.001s 74 GcV • ( A22 ) 



The time it takes to reach this flux level is 

- ~ 14 

twk 



1.7 



(A23) 



Fig. 12 shows the minimum time for GLAST to detect sources at different integral flux levels </>_8 and with different 
spectral indices from the signal and background limits, eqs. (|A13j) and (|A21j) . respectively. In order to determine the 
source's spectral index, a hardness ratio or mean photon energy (E) can be calculated for photons originating from 
a source at a known source position. For photons with energies between 100 MeV and E% — 100^2 MeV originating 
from a source with spectral index a p h, 

-1 -i a P Ji" 2 — a 

< £ > - £ '»»<?7^> ' (A24) 

Background contamination must also be considered to determine the best value and error for a p h- The upper energy 
used for u 2 depends on the nature of the source spectrum or the cutoff produced by 77 attenuation on the EBL. For low 
redshift sources, contemporaneous observations with MAGIC, HESS, or VERITAS ground-based 7-r ay telescopes could 
help c onstrain the upper energy. Higher energy analysis can enhance blazar detection. For example. iDingus & Bertsch 
(2001) reported the search for blazars in the EGRET data by looking at > 10 GcV photons. 

From the previous considerations and from Fig. 9, weak sources are best identified with GLAST by calculating 
the statistical significance S(> E)j y/B{> E) of a source as a function of photon energy E, and searching for the 
minimum in the detection significance. The actual measurement consists of the sum S(> E) and B(> E), and the 
value of B(> E) is calculated from our knowledge of the diffuse galactic and extragalactic background and GLAST 
detection characteristics. Knowledge of B(> E) will improve with our understanding of the GLAST results. For sources 
with spectra harder than the diffuse EGRB, a minimum in the combined background and signal- limited significance 
will be found by this procedure. This represents a systematic bias for GLAST to discover hard spectrum sources, and 
must also be considered when treating population statistics. At the same time, an unbiased search for sources with 
fluxes above a fixed energy E, e.g., E = 100 MeV, 1 GeV, and 10 GeV should also be made with GLAST. 
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Fig. 9. — Time required for GLAST in the scanning mode to detect a source at integral flux levels <j>—g = 1 and 0_g = 100 when 
integrating signal above energy E, for sources with photon spectral indices a p h = 2.4 (thick long-dashed curves), o p h = 2.1 (thick solid 
curves), and = 1.8 (thick dotted curves). Light dotted and dashed lines show the time required to detect a source from signal and 
background limits, respectively, for (j>—8 = 1 and a p ^ = 2.1. 

Short Timescale Flares 

GLAST reaches EGRET 2-wcck sensitivities of = 15 after ss 1 - 2 days, as can be obtained by comparing eq. 
(|A23j) with eq. (|A21j) . It reaches this sensitivity over the full sky, compared to w \/2A th of the full sky viewed with 
EGRET. For a mean photon energy of w 400 MeV measured with EGRET from a blazar with 2 < a p h < 2.5, a 
measurement of 4>s — 15 (150) corresponds to an energy flux between 100 MeV and 5 GeV of = 10~ 10 (10 -9 ) ergs 

— 9—1 

cm s . 

Consider a flare at the level of (fi-g = 200, corresponding to som e of the brightest flares seen with EGRET. Given the 
Phase 1 EGRET all-sky survey results CO ermer fc Dingusll2003[ ). we predict a flare at this level every few days, with 
large uncertainties given the small statistics of bright blazar flares. For a blazar with 7-ray spectral index a p h = 2.1, 
GLAST will have a significant 5cr detection after 



t ^ 2.6 ks 



(A25) 



using u — 1 whenever u < 1. Better sensitivity of a bright blazar flare is possible for observing times less than the 
GLAST orbital period of 1.6 hr (or 5.8 ks) if the blazar is in the FOV of GLAST. In this case, X s can be larger 
than 1/5. These timescales are of particular interest inasmuch as we expect a minimum variability timescale t var 
corresponding to a light-crossing distance 



it; 



2GM 



(1 + z) « M 8 (l + z) ks 



(A26) 



corrections should be made to this expression for a Kerr metric. Temporal power analyses of strings of data from 
bright blazars, such as 3C 279 and PKS 0528+134, could find a reduction of power at a timescale corresponding to 
the Schwarzschild radius or radius of the minimum stable orbit. Rare, bright, hard spectrum flares are of particular 
interest for GLAST analysis in order to look for 77 absorption cutoffs, whether internal (within the blob), external 
within the inner jet or galactic environment, or due to absorption on diffuse radiation fields. 

SYNCHROTRON SELF-ABSORPTION FREQUENCY 

We calculate the sync hrotron self- a bsorption (SSA) frequency in terms of the parameters of the system, l' e , T, and z. 
Using the formalism of Gould (1979), cast into dimensionless notation, we have the expression for the frequency v m 



16 



where dF v (y m ) / dv = 0, given by 

em= hj^_ = Sn_ { ^M r y eK] 2/( P+ 4) 5 ( P+2) /( P+ 4) (B1) 



and c(p) and t m are numbers of order unity (Gould 1979), otf = 1/137 is the fine structure constant, and K' e = V^k' 



From the expression for synchrotron radiation, eq. (flQ 

tI 
9 



f<"m 

C = (^)« 3 ^ 1)/2 . (B2) 



It follows from these expressions 

i/(GHz) _ (J D r 187rc(p) 



1.2 x 10 11 1 + 2; L a/i m co"tUb ct . 



(l£) 2 ] 3/(P+4) & 1/(p+4) ■ (B3) 

r 6 



We obtain 

0.22 (VEJ/ru) B(G) 1 / 6 , p = 2 



!/(GHz) = | 9.5 (/ 42 /r? 4 ) 3/7 5(G) 1 / 7 , p = 3 , (B4) 



1900 (/ 42 /r? 4 ) 3/8 i?(G) 1 /8 , p = 4 



where r' b = 10 14 ri4 cm and l' e — 10 42 Li2 ergs s _1 sr _1 . BLs are likely to be heavily self-absorbed at 5 GHz, using 
parameters So = 4, p = 3, Z42 = 1, ?T4 = 1, and B = 0.1 - 1 G. FSRQs could also be self-absorbed at 5 GHz, using 
5b — 10, p > 3, Z42 — 0.01, ri4 = 100, and B>1 G. Here the values of B are typical of equipartition magnetic field 
strengths. If GLAST shows that FSRQs vary on hour or subh-our timescales, then the synchrotron emission of blazars 
is likely to be strongly self-absorbed. 
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TABLE 1 

Sample of High-Confidence Gamma-Ray Blazars used in Analysis 



Catalog Name 




&4>% 


Redshift z 


Other Name b 


Classification 




3EG J0204+1458 


52.8 


26.4 


0.405 


PKS 0202+14 




3EG J0210-5055 


134.1 


24.9 


1.003 


PKS 0208-512 




3EG J0222+4253 


25.3 


5.80 


0.444 


3C 66A 


B 


3EG J0237+1635 


65.1 


8.80 


0.940 


AO 0235+164 


B 


3EG J0340-0201 


177.6 


36.6 


0.852 


PKS 0336-01 




3EG J0412-1853 


49.5 


16.1 


1.536 


PKS 0414-189 




3EG J0422-0102 


64.2 


34.2 


0.915 


PKS 0420-01 




3EG J0442-0033 


85.9 


12.0 


0.844 


NRAO 190 




3EG J0450+1105 


109.5 


19.4 


1.207 


PKS 0446+112 




3EG J0456-2338 


14.7 


4.20 


1.009 


PKS 0454-234 




3EG J0458-4635 


22.8 


7.40 


0.8580 


PKS 0454-46 




3EG J0459+0544 


34.0 


18.0 


1.106 


PKS 0459+060 




3EG J0500-0159 


68.2 


41.3 


2.286 


PKS 0458-02 




3EG J0530+1323 


351.4 


36.8 


2.060 


PKS 0528+134 




3EG J0540-4402 


91.1 


14.6 


0.894 


PKS 0537-441 


B 


3EG J0721+7120 


45.7 


11.1 


0.30 


S5 0716+714 


B 


3EG J0737+1721 


29.3 


9.90 


0.424 


PKS 0735+178 


B 


3EG J0743+5447 


42.1 


8.30 


0.723 


RX J0742.6+5444 




3EG J0828+0508 


35.5 


16.3 


0.180 


PKS 0829+046 


B 


3EG J0829+2413 


111.0 


60.1 


0.939 


OJ 248 




3EG J0845+7049 


33.4 


9.00 


2.172 


4C 71.07 




3EG J0852-1216 


44.4 


11.6 


0.566 


PMN J0850-1213 




3EG J0853+1941 


15.8 


6.90 


0.306 


OJ 287 


B 


3EG J0952+5501 


47.2 


15.5 


0.901 


OK 591 




3EG J0958+6533 


18.0 


9.40 


0.368 


S4 0954+65 


B 


3EG J1104+3809 


27.1 


6.90 


0.031 


Mkn 421 


B 


3EG J1200+2847 


163.2 


40.7 


0.729 


TON 0599 




3EG J1222+2841 


53.6 


14.1 


0.102 


W Comae 


B 


3EG J1224+2118 


48.1 


15.3 


0.435 


PG 1222+216 




3EG J1229+0210 


48.3 


11.3 


0.158 


3C 273 




3EG J1230-0247 


15.5 


4.10 


1.045 


PKS 1229-02 




3EG J1246-0651 


44.1 


29.6 


1.286 


PKS 1243-072 




3EG J1255-0549 


267.3 


10.7 


0.538 


3C 279 




3EG J1329+1708 


33.1 


19.3 


2.084 


OP 151 




3EG J1339-1419 


20.2 


11.6 


0.539 


PKS 1335-127 




3EG J1409-0745 


128.4 


23.4 


1.494 


PKS 1406-076 




3EG J1429-4217 


55.3 


16.3 


1.522 


PKS 1424-41 




3EG J1512-0849 


49.4 


18.3 


0.361 


PKS 1510-08 




3EG J1605+1553 


42.0 


12.3 


0.357 


4C 15.54 


B 


3EG J1608+1055 


62.4 


13.0 


1.226 


OS 111 




3EG J1614+3424 


68.9 


15.3 


1.401 


OS 319 




3EG J1625-2955 


258.9 


15.3 


0.815 


PKS 1622-29 




3EG J1626-2519 


82.5 


35.0 


0.786 


PKS 1622-253 




3EG J1635+3813 


107.5 


9.60 


1.814 


4C 38.41 




3EG J1727+0429 


30.2 


18.8 


0.296 


PKS 1725+044 




3EG J1733-1313 


104.8 


34.7 


0.902 


PKS 1730-13 




3EG J1738+5203 


44.9 


26.9 


1.375 


OT 566 




3EG J1744-0310 


48.7 


19.6 


1.054 


PKS 1741-03 




3EG J1935-4022 


93.9 


31.4 


0.966 


PKS 1933-400 




3EG J1937-1529 


55.0 


18.6 


1.657 


PKS 1936-15 




3EG J2025-0744 


74.5 


13.4 


1.388 


PKS 2023-07 




3EG J2036+1132 


35.9 


15.0 


0.601 


TXS 2032+117 


B 


3EG J2055-4716 


35.0 


20.9 


1.489 


PKS 2052-47 




3EG J2158-3023 


30.4 


7.70 


0.116 


PKS 2155-304 


B 


3EG J2202+4217 


39.9 


11.6 


0.069 


BL Lac 


B 


3EG J2232+1147 


51.6 


15.0 


1.037 


CTA 102 




3EG J2254+1601 


116.1 


18.4 


0.859 


3C 454.3 




3EG J2321-0328 


38.2 


10.1 


1.411 


PKS 2320-035 




3EG J2358+4604 


42.8 


20.3 


1.992 


OZ 486 




3EG J2359+2041 


26.3 


9.00 


1.066 


OZ 193 





a 4> P -s- peak flux (E > 100 McV) in units of 10~ 8 photons cm" 2 s" 1 . 13 Survey catalogs 
include PKS: Parkcs, O(A-Z): Ohio, 3C/4C: Cambridge, and TXS: Texas. c B: BL Lac 
object; no entry: FSRQ. 
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